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Abstract 

This  paper  presents  generic  models  for  the  effect  of  a  chemical 
toxin  on  cells  forming  the  tissue  of  an  organ.  The  models  are 
illustrative,  not  specific  to  organ  or  toxin.  Interactive  response  of 
within-tissue  toxin  and  cells  is  modeled:  cell  capability  to  modify 
(metabolize  or  bind)  toxin  is  represented,  as  is  the  alteration  of  that 
capability  by  toxin  presence.  Both  processes  are  represented  in  the 
context  of  simple  versions  of  the  cell  cycle. 

Toxin-cell  interaction  is  explicitly  represented  in  terms  of  inter¬ 
cell  chemical  signaling  that  encourages  replacement  or  repair  of 
cells.  The  impact  of  cell  death  by  apoptosis  and  necrosis  is 
represented  in  terms  of  within-tissue  toxin  concentration.  The 
model  is  explicitly  stochastic,  representing  inter-cell  and  toxin  input 
and  within-organ  concentration  in  terms  of  diffusion 
approximations.  Explicit  mathematical  discussion  is  given  of  dose- 
response  function  behavior  at  low  doses. 


1.  Introduction 


Application  of  laboratory  toxicology  data  to  environmental  and  human 
problems  of  risk  assessment  almost  always  requires  extrapolation  of  the  data 
from  the  experimentally-used  dose  regimen  to  the  exposure  conditions  of 
practical  concern,  and  from  the  animal  species  tested  to  the  species  of  concern 
(usually  man).  This  extrapolation  and  estimation  process  is  known  as  chemical 
risk  assessment.  The  risk  assessment  process  has  undergone  considerable 
evolution,  moving  from  a  qualitative  basis  for  decision  making  to  an  increasingly 
quantitative  basis,  and  from  the  use  of  default  assumptions  to  the  application  of 
mechanistic  mathematical  models  as  tools  upon  which  to  base  decisions.  In  the 
context  of  determining  safe  human  exposure  limits  to  potentially  toxic  chemicals, 
there  are  two  sub-tasks  to  be  accomplished:  estimation  of  low-dose  risk  in 


animals,  and  the  subsequent  conversion  of  animal  risk  estimates  to  human  risk 


estimates.  An  authoritative  survey  of  many  of  the  statistical  issues  and 
opportunities  is  given  by  Krewski  and  Franklin  (1991). 

Extrapolation  to  low-dose  effect,  the  first  sub-task,  can  be  accomplished  either 
by  assuming  that  biological  response  varies  in  a  specified  mathematical  manner 
(e.g.  probit,  logit)  with  organ  host  exposure,  or  by  using  physiologically-based 
pharmacokinetic  (PBPK)  models  to  relate  organ  dose  to  host  exposure;  cf. 
Krewski  and  Franklin  (1991),  chapters  8  and  9;  M.  Andersen,  H.  Clewell  Jr.,  and 
C.  Travis  have  been  prominent  in  this  research  area.  In  the  latter  approach,  multi¬ 
compartment  physiological  models  are  formulated  using  actual  tissue  volumes 
from  the  experimental  and  target  species  and  actual  perfusion  rates  to  provide  y01T 
for  chemical  transport  between  the  compartments.  Thus  the  pharmacokinetics  of  rJ 
high  to  low-dose  extrapolation  become  amenable  to  calculation,  and  external ;  loa 
measures  of  external  dose  or  exposure  can  be  translated  to  concentrations  in  the 
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target  organ  (internal  dose).  Target-organ  chemical  concentration  may  be 
translated  into  estimates  of  risk  if  a  suitable  biologically  acceptable  mathematical 
model  can  be  used  to  relate  chemical  presence  in  the  organ  to  harmful  outcome. 
Such  models,  called  pharmacodynamic,  exist  for  carcinogens,  in  the  form  of  the 
widely-used  linearized  multistage  approximation  for  cancer  dose/ response,  and 
multistage  clonal  expansion  models;  cf.  Moolgavkar  (1988).  These  models  are 
based  on  the  general  concept  that  chemical  alteration  of  the  cellular  genes  may 
give  rise  to  permanent,  heritable  changes  in  the  genetic  information  stored  in  the 
cell  nucleus,  and  lead  to  phenotypic  changes  in  the  altered  cells  that  ultimately 
cause  the  formation  of  malignant  tumors.  Simulation  is  sometimes  used  to  study 
stochastic  models  of  carcinogenesis  for  large  numbers  of  cells;  (cf.  Bois  (1992)). 

Unfortunately,  analogous  dose /response  models  are  not  widely  available  for 
toxic  responses  other  than  carcinogenesis.  The  mathematical  expression  of  even 
the  relatively  simple  concept  of  genetic  alteration  leading  to  cancer  involves 
significant  simplification  of  biological  reality,  and  significant  mathematical 
complexity.  This  state  of  affairs  is  exacerbated  when  one  attempts  to  describe  the 
interactions  leading  to  loss  of  cell  function  and  cell  death  in  tissues  of  a  whole 
animal.  The  multi-layered  control  and  response  systems  present  in  an  intact 
living  animal  are  poorly  understood  and  thus  have  not  yet  been  adequately 
modeled.  Nonetheless,  these  control  mechanisms  defend  against  the  majority  of 
toxic  effects  observed  as  a  result  of  chemical  exposure,  and  their  failure 
represents  much  of  what  is  observed  as  the  expression  of  toxicity. 

This  paper  describes  initial  formulations  of  models  for  cellular  response  to 
toxic  insult.  Beginning  with  a  simple  Markov  model,  successive  models  increase 
in  biological  resolution  or  refinement  as  well  as  mathematical  complexity.  Our 
model  formulation  takes  into  account  broad  aspects  of  the  current  state  of 
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knowledge  concerning  the  control  and  regulation  of  cellular  properties  in  tissue. 
The  design  of  these  models  is  not  to  provide  a  detailed  description  of  the 
response  of  specific  organs  to  toxic  insult,  but  rather  to  begin  the  process  of 
mathematically  describing  recent  significant  advances  in  knowledge  that  have 
been  made  in  cell  biology,  hormonal  regulation,  and  hormetic  control  mech¬ 
anisms  that  regulate  cell  response  to  toxicants.  Some  modifications  of  the  model 
formulations  may  be  needed  to  apply  them  to  specific  organs;  for  one  thing  the 
complex  geometry  of  organ  architecture  is  not  represented.  To  have  included 
complex,  three-dimensional  relationships  between  different  cell  types  would 
have  increased  the  mathematical  complexity  considerably.  But  even  in  its  present 
rudimentary  state  the  model  is  useful  for  stimulating  discussion  and  suggesting 
experimentation,  with  the  goal  of  furthering  the  understanding  of  cellular 
response  to  toxic  substances  including  the  effects  of  tissue-mediated  responses 
(e.g.  mitogenic  stimulation).  The  current  formulation  of  the  model,  or  extensions 
thereof,  is  suitable  for  experimental  evaluation  in  cell  culture  if  the  experimental 
conditions  are  properly  chosen  and  chemical-specific  details  are  added. 

We  provide  the  following  abbreviated  biological  background.  The  tissue 
making  up  organs  almost  always  involves  a  complex  geometrical  juxtaposition  of 
several  cell  types  having  specific  and  often  overlapping  functions.  This  tissue 
architecture  may  be  maintained  by  the  presence  of  a  nonliving  matrix  of  proteins 
collectively  known  as  a  basement  membrane.  Interactions  between  this  membrane 
and  the  cells  are  known  to  be  critical  to  its  stability  as  a  mature,  functioning 
entity.  The  whole  of  the  tissue  is  permeated  by  branching  blood  vessels,  each 
generation  of  which  is  successively  smaller;  these  serve  to  provide  a  constant, 
stable  milieu  in  which  the  cells  exist.  Nutrients,  control  signals  in  the  form  of 
specific  biomolecules,  and  xenobiotic  chemicals  are  brought  to  the  immediate 
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surroundings  of  the  cell  by  the  vasculature  in  the  tissue  and  cellular  metabolic 
products;  the  products  of  energy  metabolism  and  cellularly  derived  control 
biomolecules  are  removed  from  the  cellular  microenvironment  by  the  same 
means. 

The  state  of  the  cell  at  any  time  is  a  reflection  of  its  age,  the  summation  of  the 
control  chemicals  reaching  and  leaving  the  cell,  the  effects  of  xenobiotic 
chemicals  present,  if  any,  and  the  state  of  the  cells  surrounding  it.  Cellular 
contact  with  the  surrounding  cells  and  basement  membrane  act  to  supply 
chemical  signals  that  modulate  its  activity.  A  given  cell  may  be  (i)  nearly 
quiescent,  (ii)  active  biochemically,  i.e.,  producing  metabolites  of  absorbed 
materials  for  its  own  use  or  for  export,  (iii)  in  a  state  of  stress  due  to  shortage  or 
excess  of  biochemical  molecules,  (iv)  in  a  process  of  programmed  cell  death 
(apoptosis),  (v)  in  the  process  of  dying  from  chemical  insult  (necrosis),  or  (vi) 
dividing  to  form  progenitor  cells  in  response  to  a  need  to  replace  cells  already 
lost;  these  conditions  need  not  be  mutually  exclusive  at  any  point  in  time.  Some 
specific  cells  may  alter  or  completely  change  their  observable  characteristics  in 
response  to  chemical  signals.  The  most  well-known  of  these  cells  are  the 
pluripotent  stem  cells  of  the  hematopoietic  system. 

The  models  to  be  presented  here  represent  some,  but  by  no  means  all,  of  the 
features  mentioned.  They  provide  a  basis  for  proceeding  further. 

2.  Model  Structure 

In  the  following  sections  of  this  paper  we  present  a  sequence  of  mathematical 
models  that  represent  interactions  between  the  cells  in  an  organ's  tissue  and  a 
toxic  chemical  or  agent  (hereafter  called  toxin)  originating  externally  and  entering 
the  organ  and  coming  into  contact  with  the  cells,  possibly  in  modified  form.  The 
initial  models  postulate  a  mature  organ  with  capacity  to  contain,  at  most,  a  fixed 
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number  of  cells.  Broadly  speaking,  the  organ's  cells  attempt  to  remove  (e.g. 
metabolize  or  bind)  the  toxin,  but  are,  to  a  dose-dependent  degree,  also  affected 
by  the  toxin  concentration.  We  distinguish  between  cells  that  are  in  one  of  two 
stages  of  the  cell  cycle:  these  we  call,  oversimplifying  the  true  cell  cycle, 
functional,  and  dividing  (S-phase  and  mitotic).  Functional  cells  are  assumed  to  be 
capable  of  removing  toxin,  but  are  also  susceptible  to  premature  death  because  of 
the  toxin's  action.  Dividing  cells  are  in  the  process  of  DNA  replication  and  actual 
mitosis,  here  splitting  into  two  daughter  cells;  while  in  this  stage  the  cells  are 
assumed  incapable  of  removing  toxin.  However,  toxin  presence  is  assumed  to 
affect  the  cell  cycle,  e.g.  by  shortening  the  functional,  and  lengthening  the 
dividing,  stages.  The  model  proposed  emphasizes  inter-cellular  signaling:  when 
a  cell  dies,  i.e.  ceases  to  function,  it  effectively  requests  one  other  functioning  cell 
to  divide  to  produce  its  replacement.  It  is  biologically  plausible  that  signaling 
occurs  only  to  neighboring  cells,  but  the  present  model,  being  without  spatial 
characteristics,  does  not  literally  recognize  this  restriction.  Also,  signals  are 
assumed  to  be  specific:  one  signal  is  emitted  per  cell  death,  so  eventually  one 
replacement  occurs.  Alternative  formulations  are  sometimes  plausible  and  are 
easily  studied. 

3.  Single-Stage  Markov  Model 

We  now  present  several  alternative  simple  Markov  models  for  cell-toxin 
interaction,  emphasizing  the  representation  of  signaling.  Although  cells  can  age 
and  exist  in  different  stages  of  functionality  we  omit  an  account  of  this  at  first, 
but  return  to  it  later. 

Let  D(0  represent  the  number  of  functional  cells  in  the  organ  at  time  t,  and 
M(f)  the  number  of  dividing  (mitotic)  cells.  If  a  functioning  cell  that  dies 
immediately  induces  another  to  begin  division  (instant  signaling)  then 
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Co  =  D(0  +  2M(f)  (3.1) 

where  Co  is  the  maximum  organ  size,  here  assumed  constant;  realistic 
elaborations  are  possible.  Note  that  organ  size  is,  strictly  speaking,  D(t)  +  M(f), 
which  is  the  number  of  cells  present  at  time  /;  Co  acknowledges  that  space  for  all 
cells  must  leave  room  for  each  cell  currently  in  mitosis  to  eventually  become  two 
differentiated  cells.  We  describe  {0(0,  M(0)  as  a  birth  and  death  process  in 
continuous  time,  with  D(t)e(0, 1,  2, ...,  Co).  Typically  Co  is  large,  i.e.  of  order 
10lO_lOl2  for  mature  humans,  and  this  will  be  exploited  to  carry  out  an 

asymptotic  analysis. 

Let  T(f)  denote  the  concentration  of  toxin  in  the  organ  at  time  t.  We  represent 
{7X0}  as  a  diffusion  process  whose  drift  and  diffusion  coefficients  are  both 
influenced  by  D(0  so  as  to  represent  metabolic  or  binding  action,  hence  toxin 
removal,  and  by  the  toxin  input  to  the  organ,  t(0,  as  well  as  the  instantaneous 
concentration,  T(t),  itself. 

In  the  following  we  describe  specific  models,  beginning  with  the  most  simple 
and  (potentially)  least  realistic,  but  also  the  most  parsimonious.  We  explicitly 
include  the  effect  of  inter-cellular  signaling. 

3.1  Cellular  Signaling 

Suppose  a  cell  dies  from  some  cause:  apoptosis  or  necrosis.  The  effect  of  its 
signaling  to  other  cells  for  replacement  can  be  captured  by  postulating  that  a 
ghost  (release  of  growth  factor)  of  the  newly-dead  cell  now  exists,  the  function  of 
which  is  to  find  and  interact  with  a  neighboring  live  cell,  inducing  that  cell  to 
divide.  The  ghost  then  disappears  after  its  purpose  is  served.  One  can  thus 
represent  the  delay  in  mitotic  response  to  cell  death  by  adjusting  the  ghost's 
search  rate  for  functional  cells. 
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The  following  describes  a  possible  Markov  generator  that  includes  the 
signaling  (ghost)  influence.  Let  Git )  denote  the  number  of  signals  present  in  the 
organism  at  time  t,  and  introduce  the  constraint 

C0  =  D(t)  +  G(t)  +  2M(t).  (3.2) 

The  following  then  describes  the  possible  transitions  allowed  to  occur  in  a  small 
time  interval,  it,  t  +  dt). 


Signaling  Model 
(Markov  Generator) 


t 


t  +  dt 


Probability 


D(t),G(t),M(t ); 

m 


D(t)  - 1,  G(t)  +  W(0 
T(t)  +  S*dt  +  odWit) 
(cell  death,  signal  creation) 


X*Dit)dt 


—> 


D(f)-l,G(f)-l,M(f)  +  l 
T(t)  +  S*dt  +  odW(t) 

(cell  signaled  to  divide,  signal 
disappears) 


e*G(.t)(D(t)/Co)dt 

(3.3) 


D(f)  +  2,  G(t),  M(t)  - 1  dt 

Tit )  +  8'dt  +  orfW(0 
(two  newly-divided  cells  emerge) 


All  rates  of  cell  transition,  namely  X*,  & ,  and  fi  are  presumed  to  be  influenced 
by  current  toxin  concentration,  1{t).  In  turn,  the  latter  is  influenced  by  the 
metabolic  or  binding  capability  of  the  cells  to  remove  toxin;  the  net  mean 
increase  of  toxin  in  the  short  run  is  5*(IX0,D(f),t(f)).  Our  subsequent  analysis 
does  not,  in  principle,  depend  upon  specific  functional  forms  for  any  of  the  above 
parameters;  adequate  smoothness  and  differentiability  is  assumed  where 
necessary.  When  specific  functions  are  required  we  shall  use  these: 
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(3.4, a) 


A*(r(f))=Vl'T(') 

/(r(())=/v_"'T(')  (3.4,b) 

e'(T(t))  =  eoe~eTI-t);  >0,  (3.4, c) 

Note  that  it  has  been  assumed  that  increased  toxin  level  increases  cell  death  rate, 
decreases  the  rate  of  completing  mitosis,  and  slows  the  rate  at  which  ambient 
signals  reach  their  destinations.  The  net  effect  is  to  reduce  the  number  of 
functional  cells.  All  of  the  above  dependencies  are  hypothetical  and  illustrative 
only.  In  particular 

«*(r(l),D(/),T  I (f))  =  jl  W-V^)  (3.5) 

where  in  the  latter  T*(t)  represents  toxin  input  rate  to  the  organ,  and  the 
remaining  term  v*(T(f)D(f)/C0)/[l  +  K:*r(f)]  represents  the  toxin  concentration 

rate-limited  (Michaelis-Menton)  reduction  by  functional  cells. 

It  has  been  observed  that  the  realistic  effect  of  toxin  presence  on  A*,  ji,  etc. 

at  time  t  may  involve  the  entire  past  history  of  such  exposure,  e.g. 

A*(f)  =  £ [T(x),a<x<t);  the  explicit  function  may  be  illustrated  by  the  form 

A ,*(f)  =  exp  K  A*(T(x)dx)J.  A  further  model  enhancement  would  recognize  that 

toxin  presence  in  the  organ  may  result  in  cells  that  complete  completing  mitosis 

giving  rise  to  damaged  cells,  e.g.  possessing  DNA  adducts.  These  cells  now 

become  susceptible  to  repair. 
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3.2  Differential  Equations  For  Deterministic  Approximations 

Let  Co  denote  the  number  of  spaces/holes  that  cells  may  occupy.  Then 


C0  =  D(t)  +  G(t)  +  2M(f)  (3.6) 

since  the  cells  in  division  potentially  require  two  holes.  Assume  that  if  Co  »  1 
then 

D(t)/C0  exit),  G(t)/Co  y(0,  T(0/C0  -»  m  (3.7) 

in  probability  as  Co  — »  00  where  the  latter  functions  are  0(1).  The  following 
differential  equations  result  from  direct  manipulation  of  (3.3): 

Ml  =  _Aa(()  -  Ba(t  )y(t )  +  2^(1 -  a{t) -  y(t  ))i ;  (3.8) 

at  x. 

the  last  term  occurs  by  virtue  of  (3.6).  Also, 

^  =  (3.9) 

dt 

(3.10) 


We  re-emphasize  that  all  parameters  are  implicit  functions  of  (current)  toxin 
level,  T(f),  and  possibly  also  absolute  time,  t,  which  can  be  viewed  as  the  age  of  a 
mature  organ,  rather  than  an  individual  cell.  This  latter  step  is  not  taken  in  this 
paper. 

Note  that  in  (3.8),  (3.9),  (3.10)  the  original  parameters  of  (3.3)  must  be  scaled: 
r(T(D)  is  replaced  by  A(/J(0),  0*0X0)  by  0(#O),  ^*(T(0)  by  /x(0(O),  v*(T(0)  by 
v(/3(0),  t*  by  tCo,  and  k  by  kCo- 

In  Section  6  such  scaling  is  exploited  more  extensively  to  deduce  stochastic 
behavior  of  system  state  values. 
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3.3  Solution  Without  Toxin  Input 

Suppose  no  toxin  input  exists,  so  T(t)  =  0,  Vf  >  0.  Then  a  steady-state  solution 
to  (3.8)  and  (3.9)  may  occur:  set  derivatives  =  0  and  solve  the  resulting  equations 
to  find 


y  =  — ,  provided  0<-^<l 

0  8 

(3.11) 

^i-x/e) 

(3.12) 

2A +fi 

The  latter  implies  (i)  that  the  larger  the  signal  activity  rate,  6,  the  smaller  the 
ambient  signal  population,  and  the  more  quickly  does  signaling  induce  another 
cell  to  begin  division;  furthermore  (ii)  if  6  »  A,  the  death  rate,  signifying  highly 
efficient  signaling,  then  a  -4  p/(2  A  +  p),  which  is  equivalent  to  thinking  of  cells  as 
behaving  in  pairs:  once  division  completes  the  two  daughters  effectively  compete 
to  die  (rate  2A). 

Although  the  equations  (3.8)  and  (3.9)  are  non-linear  and  apparently  cannot 
be  solved  explicitly  one  can  obtain  an  approximate  time-dependent  solution  as 
follows: 

(a)  assume  quickly  reaches  steady  state,  so 

y(t)  =  A/0;  0  <t  (3.13) 

This  is  a  standard  and  useful  approximation  technique  often  invoked  in 
biochemical  kinetics  problems  that  is  variously  abbreviated  the  quasi-steady-state 
assumption  (QSSA)  or  the  pseudo-steady-state  hypothesis  (PSSH);  see  Segel  and 
Slemrod  (1989)  for  a  careful  exposition  of  its  rationale.  The  QSSA  approach  gives 
rise  to  the  classical  Michaelis-Menton  used  widely  in  pharmaceutical  kinetics  and 
compartment  models;  it  already  appears  implicitly  in  our  (3.5)  and  (3.10). 
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(b)  introduce  (3.13)  into  (3.8)  to  find  the  approximation  a(t) : 


da(t)/dt  =  -(2A  +  n)a  +  n{l  -  A/0), 

the  solution  of  which  is 

«*(<)= a(0)<f  +  [n{l  -  ye)/(2k  +  ) 

-» /x(l  -  A/0)/(2A  +  p)  asf->°°, 

as  in  (3.12). 


(3.14) 


(3.15) 


3.4  Dose-Response  for  Small  Steady  Toxin  Exposure 

Of  interest  in  risk  analysis  is  the  behavior  of  the  dose  response  curve  for  small 

values  of  (toxic)  dose.  We  approximate  this  by  finding  expressions  for 

da{r)  dy{r)  and  dp(r) 
dr  r  =  O'  dr  t  =  0'  dr  r  =  0 

Suppose  rit)  =  r  =  r  ~  0,  a  constant,  and  that  exposure  has  proceeded  for  some 
time  so  that  a  steady-state  condition  has  been  reached;  this  is  modeled  by  letting 

da{t)  _  dy{t)  _  d/?  _  Q  .r  (3  8)  _  (3.10).  To  study  the  dose-response  for  t=0, 
dt  dt  dt 

differentiate  (3.10)  with  respect  to  fusing  the  function  8  in  (3.5)  (drop  subscript 
"o"  for  convenience) 

dv  dp  ap  vda  p  va  dp  [  vap  dp  (316) 

0_  dpdrl  +  Kp  dr  l+Kp  l  +  tcp  dr  (i  +  Kp)2  dr’ 

Since  P  =  0  for  r  =  0,  equation  (3.6)  yields  for  the  (scaled)  rate  of  toxin 

concentration  at  low  exposure 


dr  t  =  0  var  =  0  0 


(3.17) 


This  shows  explicitly  how  low-dose  toxin  level  in  the  organ  tissue  increases  as  a 
function  of  cell-cycle  parameters  A  and  \i,  signaling  efficiency,  0,  and  toxin 
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metabolic  or  binding  rate  v.  All  of  these  parameters  are  evaluated  at  very  low 
(zero)  toxin  levels. 

Next  consider  steady-state  rate  of  change  of  scaled  functional  cell  fraction, 

—  =  aT  =(Xr  pr;  this  is  the  slope  of  a  dose-(functioning  cell  count)  response 

dx  p 

curve  at  very  low  levels  of  toxin.  If  (3.12)  is  differentiated  and  some 
rearrangements  made,  and  if  y=  A/0  as  in  (3.11),  then 

=  [-%^(2 + A/^  +  2(1  -  r }) + (2A(1  -  y)) + (2A  +  >x)(pi/e)y]  •  (318) 

Since  (3r  >  0,  the  sign  of  oy  is  determined  by  the  above  linear  combination  of  the 
derivatives  of  the  cell  cycle  and  signaling  parameters  A,  /i,  and  0  at  T,  and  hence 
p,  equal  to  zero.  Note  that  if  the  illustrative  relations  (3.4)  are  invoked  as  given, 
then  ap  <  0,  hence  <  0  so  a  small  increase  in  toxin  induces  a  decline  in  the 
number  of  functioning  cells;  this  is  intuitively  acceptable.  However,  a  change  in 
sign  of  any  of  the  parameters'  derivatives  is  capable  of  reversing  the  sign  of  the 
slope  of  the  dose  response  curve,  thus  representing  a  tendency  for  hormesis  at 
low  dose. 

4.  Models  with  Cell  Aging 

In  this  section  we  generalize  the  models  of  the  last  section  to  allow  cells  to 
age,  and  hence  exist  in  possibly  different  stages  of  functional  effectiveness. 

We  restrict  attention  to  just  two  age  stages,  which  we  term  new  (n)  and  old  (o), 
distinguishing  them  by  subscript.  An  arbitrary  number  of  such  age  states  can  be 
introduced,  but  at  the  cost  of  increasing  the  number  of  parameters. 

Let  Dn(t)  (respectively  Do(0)  represent  the  number  of  new  (respectively  old) 
cells  at  time  f.  M(t)  represents  the  number  of  dividing  (mitotic)  cells  at  time  t  and 
G(t)  the  number  of  ambient  signals  or  ghosts  as  before.  The  spatial  constraint  is 
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(4.1) 


C0  =  D„(0  +  D0it)  +  Git )  +  2M(f). 
Here  are  the  transition  probabilities  assumed: 

Signaling  Model  with  Age  Dependence 
(Markov  Generator) 


t 

t  +  dt 

Probability 

Dnit),  Do(t),  — » 

Dnit)  - 1,  Do(0,  G(f)  +  l,M(f) 

XnDniMt 

G(t),M(ty, 

T(0  +  5*dt  +  odWit ) 

Tit) 

(new  cell  death,  signal  creation) 

Dn(f),  D0it)  - 1,  G(f)  +  l,M(f) 

X0D0it)dt 

Tit)  +  8*dt  +  odWit) 

(old  cell  death,  signal  creation)! 

—> 

D„(f)  - 1, 00(f),  G(0  -  l>f(0  +  1 
T(0  +  5 *dt  +  odWit) 

(new  cell  signaled  to  divide, 
signal  disappears) 

6*nGit)iDnit)/Co)dt 

-> 

Dnit),  D0it)  - 1,  Git)  - 1,  Mit)  +  1 
Tit)  +  S*dt  +  odWit) 

(old  cell  signaled  to  divide, 
signal  disappears) 

eoG(t)(Do(t)/C0)dt 

-» 

Dnit) +  2,  D0it),  Git),  Mit) -l 

Tit)  +  8*dt  +  odWit) 

(two  newly-divided  new  cells 
emerge) 

tfMit)  dt 

-» 

D„(0-1,Do(0  +  1,  Git),  Mit) 
Tit)  +  8'dt  +  odWit) 

(a  new  cell  becomes  an  old  cell) 

<(>  *Dnit)  dt 

Once  again  all  rates  of  cell  transition,  namely  Aj,,  A*,  6*n,  0*,  ji  are  presumed 
to  be  influenced  by  current  toxin  concentration.  In  turn,  the  latter  is  influenced  by 
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the  metabolic  capability  of  the  cells  to  remove  toxin,  S*(T(t),  Dn(t),  D0(t), 
When  specific  functions  are  required,  we  shall  use  for  illustration 

S'(T(t),Dn(l),DJt),  z(t))=z  (l)-v'  -  vl  T('f 1 (4.3a) 

l  +  Kfllyt)  l  +  Kplyt) 


^  J  ) 


p  +  q  =  1,  p^O  (4.3b) 


(4.3c) 


(4.3d) 


(4.3e) 


The  form  of  (4.3b)  permits  the  adjustment  of  organ  toxin  concentration-caused 
cell  death  rate  to  be  flexibly  adjusted;  if  £n\  <  0  the  inner  bracket  must  have  its 

sign  reversed:  (  )+.  As  before,  the  above  forms  are  hypothetical. 

4.1  Differential  Equations  for  Deterministic  Approximations 

Let  Co  denote  the  number  of  spaces/holes  that  cells  may  occupy.  Then 


C0  =  Dn(t)  +  Do(t)  +  G(t)  +  2M(t). 


Assume  that 


->  «n(0/  Y(t)>  ^  P(t)  (4.4) 

Co  C0  Cq  C0 

in  probability  as  Co  -» °°.  The  following  differential  equations  result  from  direct 

manipulation  of  (4.2)  with  scaled  rate  functions 

= -\k + <t>]an(t)  -  @nan  (0/(0 + 41  -  «»(0  -  ao  (0  -  y(0]  <4-5) 


<fro(0  _ 


V*o(O  +  0«n(O  ^o®o(0/(0 
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~  -[8nan(0  +  eo“o(0]r(<)  +  +  A<j“o(') 


(4.7) 


(rt-vdll-® _ yg(t)  M  (- 

■if  W  V"“"l,l  +  ^(f)  Vo“oWl  +  ^(f) 

These  equations  can  be  solved  numerically,  but  not  in  closed  parametric  form. 


4.2  Solution  Without  Toxin  Input 

Suppose  no  toxin  input  exists  so  T(t)  =  0,  t  £  0.  Then  a  steady-state  solution  to 
(4.5)  -  (4.7)  may  occur;  set  derivatives  equal  to  0  obtaining  the  following 
equations. 

0  =  -(An  +  <p)an  -  enanY  +  n[l  -an-a0-y]  (4.9) 


0  =  -X0a0  +  <t>an-eoa0Y  (4.10) 

0  =  [— —  0o«o]y  +  +  A0tt0.  (4.11) 

Equation  (4.10)  yields 


=  (l/b)[Ao  +  0oy]ao.  (4.12) 

Substituting  (4.12)  into  equation  (4.11)  yields 

0 = HMVtfpo + M  -  e0]y + +M + a„)  (4.13) 

after  dividing  through  by  Oq. 

Rewriting  (4.13)  we  obtain 


o=(-e„e0/^)y2-r[e»(V^+e0-^(V#o]+[(4/«)+iK.  (4.u> 

Thus  y  satisfies  a  quadratic  equation.  Let  Yo  be  the  positive  solution  to  the 
equation.  If  A«(0)  =  0,  as  in  (4.3b) 
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If  Yo  <  1,  then  y=  y0.  If  y0  >  1  then  ambient  signals  and  not  functional  cells 
dominate  the  organ,  which  in  reality  would  be  long  defunct  even  without  toxin. 


6n 

If  ~t  <  1  and 
<? 


e0  L  0 J/ L 


(4.16) 


then  y0  <  1. 

To  find  do  equation  (4.9)  is  used.  Putci  =  (1/ <t>)[%o  +  60$.  Then 


rfl  -  Yo]  =  {(4  +  <P)C\  +  en^i  +  rfci  +  fl}a0  =  {2A0  +  //[q  + 1]  +  2Vl}«o- 


Thus, 


ao  =  rf1  ”  7o]/{22v  +  /41  +  cl]  +  2Vl} 


(4.17) 


«n  =cia0. 


(4.18) 


4.3  Dose-Response  For  Small  Steady  Toxin  Input 

As  mentioned  earlier,  the  behavior  of  an  organ-level  dose-response  curve  for 
small  values  of  toxic  dose  is  of  interest  in  risk  analysis.  In  this  section  we  indicate 
that  the  derivative  of  various  cellular-level  responses  with  respect  to  r  may  be 
evaluated  at  t=  0. 

If  t  =  0  then  (3  =  0.  Set  the  left-hand  side  of  (4.8)  equal  to  zero,  then 
differentiate  the  right  side  with  respect  to  r.  After  setting  p  =  0,  this  yields 


Thus 


0  =  1“  \ynan  +  voao]fa- 


drr  =  0  vnan  +  v0a0 


(4.19) 


(4.20) 


as  is  physically  plausible. 
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It  may  be  shown  by  further  differentiation  that  under  the  specific 
parameterization  of  (3.4),  or  suitable  alternatives,  all  derivatives  of  long-run 
cellular-level  response  can  be  explicitly  evaluated  at  t=  0,  as  was  shown  earlier. 
The  complex  details  are  omitted  here. 

5.  A  Model  with  Cell  Aging  and  Organ  Tissue  Growth 

In  this  section  we  generalize  the  previous  model  to  allow  cells  to 
spontaneously  enter  mitosis  and  multiply,  thus  effectively  increasing  the  organ 
size,  where  the  latter  is  defined  as  the  number  of  cells  in  existence,  in  either  state. 
In  the  previous  section  cell  replication  only  occurred  to  replace  other  cells  that 
previously  have  died,  possibly  as  a  result  of  toxin  action. 

Consider  the  region  to  be  occupied  by  cells  to  consist  of  spaces /holes  that 
may  be  filled  by  cells.  Let  Co  be  the  maximum  number  of  such  holes,  in  this  case 
the  number  in  an  essentially  mature  organ.  Initially  not  all  such  holes  are  active, 
for  the  organ  is  immature  and  hence  growing.  Hence  let  C(t)  be  the  number  of 
active  holes  at  time  t,  while  Co  -  C(f)  are  currently  quiescent  or  inactive.  Now, 
generalizing  the  previous  setup,  let  the  current  number  of  active  holes/ spaces  be 

C(t)  =  D„(f)  +  D0(t)  +  Git)  +  2 Mit)  (5.1) 

and  augment  the  previous  model  to  allow  C(0  to  gradually  grow,  although 
remaining  bounded  by  Co- 

5.1  Organ  Growth  Stimulated  by  Signaling 

In  this  model  the  presence  of  inactive  holes /spaces  effectively  encourages  the 
active  cells,  D„(f)  and  D0(f),  to  increase  the  number  of  signals  and  hence  the 
number  of  active  spaces  capable  of  accomodating  cells.  We  are  implicitly 
modeling  the  competing  effects  of  positive  and  negative,  or  inhibitory,  growth 
factors. 
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The  model  is  specified  by  the  following  transitions. 

Signaling-Driven  Organ  Growth 
(Markov  Generator) 


t 

t+dt 

Probability 

Dnit),  D0(t),  — ¥ 
G(t),  M(t), 

at)  nt ) 

Dn(t),  Doit),  Git)  + 1 ,  Mit),  at)  + 1 
T(0  +  Sdt  +  odWit) 

(a  quiescent  cell  space  becomes 
active,  signal  creation) 

-4 

Dn(t)  - 1,  D0(t),  G(t)  + 1,  M(0,  at) 
Tit)  +  Sdt  +  odWit) 

(new  cell  death,  signal  creation) 

^ Dnit)dt 

-4 

Dnit),  D0it)  - 1,  Git)  +  1,  Mit),  at) 
Tit)  +  Sdt  +  odWit) 

(old  cell  death,  signal  creation) 

X*0D0(t)dt 

-4 

Dnit)  +  2,  D0it),  Git)  +  1,  Mit)  - 1,  at) 
Tit)  +  Sdt  +  odWit) 

(two  newly-divided  cells  emerge) 

tfMit)dt  (5.2) 

-4 

Dnit)  - 1,  Doit)  +  1,  Git),  Mit),  at) 
Tit)  +  Sdt  +  odWit) 

(new  cell  becomes  old) 

<fDnit)dt 

-4 

Dnit)  - 1,  D0it),  Git)  - 1,  Mit)  +  1,  C(f) 
T(0  +  Sdt  +  orfW(0 
(new  cell  signaled  to  divide, 
signal  vanishes) 

<*<(! 

-4 

D„it),  D0it)  - 1,  Git)  - 1,  M(f)  +  1,  C(f) 
Tit)  +  Sdt  +  odWit) 

(old  cell  signaled  to  divide. 

»mdt 

signal  vanishes) 
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♦  ♦  ^  1^  to*  to*^ 

Once  again,  all  rates  of  cell  transition,  namely,  A„,  A0,  v,  <p ,  H,  gn,  g0  are 
presumed  to  be  influenced  by  current  toxin  concentration.  When  specific 
functions  are  required  we  shall  use 

(5.3) 

.  cm)=^o‘<m 

in  addition  to  those  of  (4.3  a-e). 


Once  again  assume  that  as  Co  -» °° 


in  probability. 


m  7^ ->»*<).  (5.4) 

Cn  Lfi 


The  following  differential  equations  result  from  direct  manipulation  of  (5.2). 

=  -[^n  +  0]a»  (0  ■ -  en  y  (0  +  Mn(0  -  «n  (0  "  «o  (0  "  /(0)  <5-5> 

^  =  -A0«0(f)  -  e0  ^y(f) + **„(*)  (5-6) 

Mi)  =  _  6n  0  °^  y(t)  +  Anan(t)  +  A0a0(f)  +  [&ao(0  +  |„a„(f)][l - *?(*)]  (5-7> 

dt  T)  (t) 

M2  =  [fca„(f)+&a„(0][l-»W]  (5-8) 

(5.9) 


6.  Stochastic  Differential  Equation  Models 

The  size  of  the  state  space  makes  the  time-dependent  behavior  of  the 
continuous-time  Markov  chain  models  of  Sections  3-5  difficult  to  study.  One 
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approach  to  studying  the  behavior  is  by  Monte  Carlo  simulation;  cf.  Bois,  et  al. 
(1992).  Another  approach,  adopted  here,  is  to  approximate  the  continuous-time 
Markov  chain  model  by  a  diffusion  process.  Note  that  the  continuous-time 
Markov  chain  model  has  absorbing  states;  for  example,  any  state  with  Dn(t)  + 
D0(t)  =  0  is  absorbing.  Of  course  no  self-respecting  organ  would  ever  start  life  in 
such  a  state,  and  a  living  organ  would  presumably  die  from  other  unmodeled 
causes  long  before  such  a  state  is  reached,  i.e.  when  D»(t)  +  D0(t)  <  D,  some  lower 
limit.  Barbour  (1976)  discusses  how  long  and  over  what  ranges  the  underlying 
continuous-time  Markov  chain  is  approximated  by  a  diffusion  process  of  the  type 
we  derive;  see  also  discussion  to  McNeil  and  Schach  (1973).  Stochastic 
differential  equation  models  can  be  written  for  all  the  models  of  Sections  3-5. 
We  will  illustrate  our  approach  by  writing  down  a  system  of  stochastic 
differential  equations  for  the  multivariate  process  (D„(f),  D0(t),  G(0,  C(t),  T(t ))  for 
Markov  generator  (5.2).  The  system  of  stochastic  differential  equations  is  as 
follows;  as  explained  subsequently,  the  dW- terms  are  normally  distributed. 

-<t>'(T(t))Dn(t)dt  +  f/ (T(t)){C(t) -  Dn(t) -  D0(t ) - G(f))dfJ\in  (f'f  +  dt )  ) 

» 

- - - - — -  ,  .  ( diffusion = random^ 

-^£9* (T(t))G{t)(Dn (f )/ C{t)) dW8n  (f)  [change  in  (t,t  +  dt) ) 

Je4  H*(T(t))(C(t)  -  Dn(t)  -  D0(t)  -  G(t))±  dWM(t) 
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dD0(t)  =  -£0(T(t))D0(t)dt  -  0„  (T(t))G(t)(D0  (t)/C(t))dt 
+0  {T(t))Dn(t)dt 


(drift) 


~^<{T(t))D0(t)  dWXo  (t)  +  p*  (T(t))Dn(t)  dW$  ( t ) 


~‘J£&o(T(t))G(t)(D0(t)/C(t))  dWffg  (t) 


(6.2) 


^diffusion) 


dG(t) = -e'„(T(t))G(tXD„(i)/c(i))di  -  e0(m)G(,iD0(,)/c(t))di 
+X0(T(t))D0(t)dt  +  X"„(T(t)p„(f)dl 


+£(T(l))D„(t) 


1- 


Q> 


+  i(T(‘))D0(t) 


1_£W 

.  Co  J 


(drift) 


-tKc mmD„(o/c(t))dwejt ) 


-^m)G(t)(D0  (l)/C(0)  dWeo  (f ) 


^eXi(T(t))D0(t)  dWK  (t) + p:(T(t))D„(t)  dWK  (t) 

rfWjo(f) 


1 


£^(r(i))D0(() 


(6.3) 


(diffusion) 
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^-  =  £(T(t))D0(t)  1-^  +£(T(t))Dn(t)  1-^  j(drift) 


(6.4) 


^  =  5,(r((),D„((),D„(f),'r‘(l))<«}(drift) 

+aT(T(t),D„(l),D0(t),r(t))dWaT(t) 

+JCoOrdWr(t) 

where  {w^  (0}/{wAo  (O}'{w0„  (O}'{w0o  (0}/{w/i(0}/{w^o(0}/{w^  (0}'{wctt  (O}/ 
and  (w^f)}  are  independent  standard  Brownian  motions;  each  dW(f)-term  is 
thus  Gaussian  with  mean  zero  and  variance  dt. 

Setting  the  constant  e  =  1  lets  the  variances  of  the  change  in  the  number  of 
active  cells,  ghosts  in  identified  states,  and  active  spaces  be  equal  to  the  mean 
respective  changes;  this  represents  Poisson  variability.  Setting  the  constant  £  >  1 
implies  that  the  variability  is  larger  than  Poisson  variability.  This  additional 
variability  may  be  the  result  of  inhomogeneities  in  the  organ  that  are  not 
explicitly  modeled;  cf.  Bass,  Robinson  and  Bracken  (1978)  for  discussion  of  a 
distributed  liver  model  with  random  inhomogeneities.  Such  overdispersion  is 
frequently  encountered  in  practice,  and  may  be  accounted  for  parsimoniously  as 
we  do;  see  McCullagh  and  Nelder  (1983). 

We  assume  that  as  Co  00 


(6.5) 

(diffusion) 


^-»o»  M,  ^ 
Co  c0 


a0{t), 

Co  Cq 
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in  probability  and  rewrite  Dn(t),  D0(t),  G(t),  C(t),  and  T(t)  as  follows. 

Dn(t)  =  CoaM+^X„(t)  (6.7a) 

D0(t)  =  C0a0(t)  +  ^X0(t)  (6.7b) 

'  G(l)  =  C0y(t)+^Xg(t)  (6.7c) 

C(t)  =  C0r,(t)  +  ^Xc(t)  (6.7d) 

T(t)  =  CoP(t)  +  ^Y(t)  (6.7e) 

dr  (l)=C0T(l)+^<rTdWI(t)  (6.70 

The  transition  rates  will  be  scaled  and  expanded  as  follows: 

4(T(f))  =  A„(T(f)/C0)  =  Um)  +  K{mW*)NCo  +  0(1)  (6.8a) 

&(m)=Km/co)=K(m)+K(m)m/-Jco <6-8w 
M'(r(f))=p(r(()/co)=M(«»))+/i'(«0)n<)/1/Q  +0(1)  (6.8c) 

e;(T(f)) = e„(T(t)/c0) = e„(/J(f))+ e'n(mM‘)N5>  +°(i)  <6-8d> 

<£  (r(<)) = e„(T(()/c0) = e„(«f))+ e;(«f))  +o(i)  (6.8e) 

£(r(t)) = €„(r(t)/c0) = !„(«())+ aWOinO/VQ  +  o(i)  (6.80 

fe*(rw) = &(t(o/q>) = &wo) + «(«<))  v«A/q + om  (6.8g> 
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(6.8h) 


S'(Dn(t),D0(t),m  T  (())  =  C08(d„(I)/C0  ,D„(f)/C0  ,T(f)/C0 ,  T  (i)/C0) 

=  C06(a„(t),cco(t),p(t),t(t)) 

«  M)+r^ i-  Ml+r^iM+om 
+c°5“”  Co5“»  Vq  °^Vq  () 

(TrfDJO-DoW.n'^'t'^^TfDnfO/Co.D.flVCo.nO/Cb-^fO/Co) 

(6.8i) 

= V^oCTrK(0/«o(0/^(0^(0) + OW 

Dividing  equations  (6.1)  -  (6.5)  by  Co,  scaling  and  substituting  (6.7a)  -  (6.7f) 
results  in  a  system  of  equations  involving  the  deterministic  components  an(t)r 
a0(t),  y(t),  rj(t),  and  /3(f)  and  the  corresponding  stochastic  components  Xn(t),  X0(t), 
Xg(t),  Xc(t),  and  Y(t).  See  Appendix  1  for  details. 

6.1  Second  Moments 

The  stochastic  differential  equations  (A.l)  -  (A.10)  of  Appendix  1,  derived  for 
the  stochastic  disturbances,  can  be  written  in  matrix-vector  form  as 

dZ(t)  =  A(f)Z(f)  +  B(f)  rfW(f)  (6.9) 

where  in  the  present  example  Z  is  a  5x1  column  vector,  A  is  a  5x5  matrix,  dW(f)  is 
10x1  column  of  independent  Gaussian  noise  terms  and  B(f)  is  a  5x10  matrix. 

It  follows  from  theorem  (8.5.5)  of  Arnold  (1974)  that 

-f  e[z ;?(()]  =  2XA,y(l)E[zi(l)Z/(f)]  +  X(B,yW)2  (6.10) 

at  .  j 

and 
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1e[z,(!)Z,(!)]  =  X  Aa(t)E[zfc  (t)Zy  (f )] 
at  k 

+XA,i(l)E(Zt(l)Zi(l)]  (6.11) 

k 

+XBil(OB,it(l) 

It 

Thus,the  approximate  distribution  of  At(f),  (respectively  Do(0/  G(t),  C(t),  T(t ))  is 
normal  with  mean  Coan(t),  (respectively  CoOo(t),  Co)(0  CoV(0,  CofKt))  and  variance 
(respectively  C0e[Xq  (f)],  C0E[x|(f)],  C0Ejxc2(f)j,  C0E[y2(t)]).  The 
approximate  covariance  of  Dn(t)  and  D0(t)  is  CoE[Xn(f)X0(0],  etc. 

We  point  out  that  exactly  equivalent  results  can  be  obtained  by  introducing 
transforms  for  the  state  variables,  computed  by  taking  conditional  mathematical 
expectations  of  exponential  functions  of  the  Markov  generators.  This  is  the  route 
followed  by  McNeil  and  Schach  (1973);  see  also  Carpenter,  Gaver,  and  Jacobs 
(1993).  The  transform  approach  provides  a  verification  that  the  asymptotics 
described  here  actually  lead  to  limiting  Gaussian  distributions  and  to  Ornstein- 
Uhlenbeck  processes.  See  Ethier  and  Kurtz  (1986)  for  a  more  rigorous  discussion. 
The  approach  taken  here,  while  heuristic,  is  intuitively  appealing  and  reaches  the 
same  ultimate  conclusion. 

6.2  Numerical  Examples 

In  this  section  we  present  three  numerical  examples. 
a.  Constant  Rate  of  Input  of  Toxin 

We  assume  chronic  dosage  of  the  organ  by  toxin  that  is  delivered  there  at 
constant  mean  rate  but  with  substantial  randomness  around  that  mean. 

Figures  1-5  depict  the  time  variation  of  the  expected  numbers  of  new  cells, 
of  old  cells,  of  total  number  of  active  (new  plus  old)  cells,  number  of  signals 
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(ghosts),  and  amount  of  toxin  plus/minus  two  population  standard  deviations. 
There  is  no  organ  growth.  The  initial  conditions  are  the  limiting  moments  with 
no  toxin  input.  The  toxin  input  rate  0.2  units/cell  for  all  time.  The  new  cells 
without  toxin  have  mean  time  100  before  transitioning  to  old.  The  old  cells  have 
mean  lifetime  of  100  without  toxin.  The  number  of  cells  in  the  organ  is  1.5x10s. 
The  organ  size  is  about  that'  of  a  mouse  liver.  Since  the  cells  in  a  mouse  liver  have 
a  mean  lifetime  of  about  200  days;  cf.  Bois,  et  al.  (1992),  the  above  is  consistent 
with  mouse  liver  behavior.  The  MATLAB  4th  and  5th  order  Runge-Kutta- 
Fehlberg  numerical  integration  method  was  used  to  obtain  the  solutions.  The 
MATLAB  plotting  algorithms  were  used  to  produce  the  figures. 

There  is  extra  Poisson  variability:  e  =  104. 

Note  that  the  total  mean  number  of  active  (new  and  old)  cells  decreases  as  the 
result  of  toxin  input,  as  seems  intuitively  reasonable.  The  mean  number  of  new 
cells  initially  decreases,  but  eventually  increases  to  a  new  steady-state  value 
which  is  larger  than  the  value  with  no  toxin.  The  mean  number  of  old  cells 
initially  increases  briefly,  but  then  decreases  to  a  new  steady-state  value  below 
the  value  for  no  toxin.  The  mean  amount  of  toxin  and  the  mean  number  of 
signals  (ghosts)  both  increase  to  new  steady-state  values. 

The  above  hypothetical  example  illustrates  how  model  parametric  inputs 
translate  into  an  account  of  the  dynamics  of  subpopulations  of  cells.  The  details 
of  the  transition  from  one  steady-state  situation  to  another  are  of  possible  interest 
in  that  the  organ  may  be  in  jeopardy  during  that  transitional  period  because  of, 
for  example,  new  cell  initial  downward  fluctuation.  The  effect  on  organ  mortality 
is,  however,  not  explicitly  modeled. 
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b.  A  Toxin  Pulse  or  Bolus  Dose 

Figures  6-10  display  the  mean  number  of  new  cells,  old  cells,  active  cells 
(new  and  old),  signals  (ghosts)  and  amount  of  toxin  plus  and  minus  2  standard 
deviations  for  the  case  in  which  the  rate  of  toxin  input  is  0.2  per  cell  for  the  first 
50  time  units  and  0  thereafter.  Once  again  there  is  no  organ  growth.  The  initial 
conditions  are  steady-state  values  when  there  is  no  toxin  input.  Note  that  the 
number  of  active  cells  decreases  during  toxin  input,  then  increases  and  slightly 
overshoots  the  steady-state  value  with  no  toxin  before  returning  to  its  steady- 
state  value  with  no  toxin. 

c.  Organ  Growth 

The  model  of  Section  5.1  and  (6.1)  -  (6.5)  is  initialized  with  the  steady-state 
values  for  the  second  moments  with  no  input  of  toxin.  However,  only  1/2  of  the 
available  spaces  are  active  (may  contain  cells);  the  steady-state  mean  of  new  and 
old  cells  and  ghosts  for  no  toxin  input  are  multiplied  by  1/2  and  used  as  initial 
values.  Figures  11-15  display  the  mean  numbers  of  new,  old,  and  active  cells, 
ghosts,  and  active  spaces  plus /minus  two  standard  deviations  for  the  case  of  no 
toxin  input  and  parameters  £ n,0  =  0.5,  £n,i  =  0.5,  £o,o  =  1,  £ o,i  =  0.5.  The  mean 
number  of  active  spaces  is  the  maximum  number  by  time  20.  However,  the  mean 
number  of  old  and  new  cells  is  still  adjusting  by  time  100.  The  mean  number  of 
old  cells  initially  decreases  before  increasing.  The  mean  number  of  new  cells 
initially  increases  before  decreasing. 

The  present  model  might  be  useful  for  describing  the  aftermath  of  a  partial 
hepatechtomy:  about  one-half  of  a  mature  liver  remains,  and  grows  back  with 
subpopulations  of  cells  responding  as  shown.  Note  that,  according  to  the  present 
model  in  Figure  12,  the  organ  must  withstand  an  early  signal  to  old  cells  to 
undergo  substantial  depletion  to  enter  mitosis  and  bring  forth  replacements;  the 
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latter  are  the  new  cells  depicted  in  Figure  11.  It  may  be  that  the  sudden  depletion 
of  old  cells  to  their  low  point  at  about  t  =  6  (Figure  12)  actually  puts  the  organ  in 
jeopardy. 

Figures  16  -  21  display  the  mean  number  of  new,  old,  and  active  cells,  ghosts 
and  active  spaces  plus/minus  two  standard  deviations  for  the  same  parameters 
as  used  in  Figures  11  - 15  but  with  two  levels  of  toxin  input.  The  dashed  lines 
correspond  to  no  toxin  input.  The  solid  lines  correspond  to  a  constant  toxin  input 
of  0.2  per  cell  per  unit  time.  The  graphs  with  positive  toxin  input  have  the  same 
general  shape  as  those  with  no  toxin  input.  The  mean  number  of  new  cells  is 
initially  less  for  a  positive  toxin  input  than  for  no  toxin  input.  However,  by  time 
100  the  mean  number  of  new  cells  with  positive  toxin  input  is  becoming  larger 
than  that  for  no  toxin  input. 

The  effect  of  positive  toxin  input  on  the  mean  number  of  old  cells  is  that  it  is 
initially  larger  than  that  for  no  toxin  but  by  time  100  the  mean  number  of  old 
cells  with  positive  toxin  is  less  than  that  for  no  toxin.  The  effect  of  positive  toxin 
input  on  the  mean  number  of  active  (old  and  new)  cells  is  to  decrease  it  and 
delay  its  approach  to  a  steady  state  value.  The  mean  number  of  ghosts  is  larger 
with  positive  toxin  input.  Positive  toxin  input  delays  the  overall  growth  in  the 
number  of  active  spaces. 
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mean  number  of  new  cells 


Figure  1 
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mean  number  of  old  cells 
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Figure  2 
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mean  number  of  active  cells  plus  minus  2  standard  deviations 
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mean  number  of  new  cells 


mean  number  of  new  cells  plus  minus  2  standard  deviations 


Figure  6 
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mean  number  of  old  cells  plus  minus  2  standard  deviations 
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active  cells 


mean  number  of  active  cells  plus  minus  2  standard  deviations 
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mean  number  of  new  cells  plus  minus  2  standard  deviations 


Figure  11 
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mean  number  of  old  cells 


Figure  12 
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mean  number  of  active  spaces 


Figure  14 
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mean  number  of  ghosts 
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mean  number  of  new  cells 
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0.5  of  the  spaces  are  active  initially 

mean  number  of  new  cells  plus  minus  2  standard  deviations 


Figure  16 
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mean  number  of  old  cells 
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0.5  of  the  spaces  are  active  initially 

mean  number  of  old  cells  plus  minus  2  standard  deviations 
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mean  number  of  old  and  new  cells 
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mean  number  of  ghosts 
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0.5  of  the  spaces  are  active  initially 

mean  number  of  ghosts  plus  minus  2  standard  deviations 
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mean  number  of  active  spaces 
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0.5  of  the  spaces  are  active  initially 

mean  number  of  active  spaces  plus  minus  2  standard  deviations 
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solid=tau=(.2/1 .5)*1 0**-8 
dash=tau=0 
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Figure  20 
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mean  amount  of  toxin 


0.5  of  the  spaces  are  active  initially 

mean  amount  of  toxin  plus  minus  2  standard  deviations 


lam  n0= 1  /1 00;  lam  o0= 1  /1 00 
mu0=25;mu1=.1 
phi0=1/100;  epsilon=10000 

thetan0=thetao0=1  ;thetan  1  =thetao1  =.5 
kapn=kapo=(1/1 .5)*10**-8 
nun=nuo=.4*(1/1 .5)*10**-8 
lamnl  =lamo1  =theta1  =phi1  =(0.5/1 .5)*1 0**-8 
sig2tau=(1 0000*tau/1 .5)*(1 0**-8) 

sig2tox=.4*1 0000*(alphan+alphao)*beta/(1  +beta) 

tau=(2/1.5)*10**-8 
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Figure  21 
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Appendix  1 

This  appendix  provides  the  detailed  asymptotic  development  of  both  the 
deterministic  approximation  (5.5)  -  (5.9)  but  also  for  the  stochastic  components 
introduced  in  (6.7).  Start  by  introducing  the  normalization  (6.7)  into  (6.1)  -  (6.5) 
and  divide  by  Co-  The  following  results  are  obtained: 

<fa„(()+-4-dx„(f)=-  um)+Km)^ I“»(')+^x-.(f)  * 

-  «>(/?(<)) +  f  (/5(f))#  a„(f)  +  -ji-X„(f)  A 
+  /i(/3(f))  +  /i'(0(f))-^5=  |[j)(f)-a„(f)-a0(f)-y(f)] 


-e„(«f))+fl;(W)#  r(0+#xs(f) 


(0-x„(()-xo(f)- 


--^L-^n(y3(f))a„(f)ilVAi,(f)--^-X«f))a„(f)dW0(l) 


(A.l) 
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dy(t)  +  lprdXg(t)  = 
VLo 


vff'iir  1  1  a”^+  A:  x«(0 

-  ^  W0)+ W))^[r(0 + ^f)+^xc(() 


v/aftW  ^(0 


«a(0  +  -7=Xo(f)' 

- 2g° -  dt 


-  »<,(«*)) +«(«0)_^  I  rW+^^wj  ,  Xi 


«nW)+6(«0)^  a.(0+^4-^)-^rxe(0  * 
umhum)^  «»(<)+ j3-'x»(')  i-fi(<)-^xc(()* 


+  XoimhKW)^  “»(f)+^=x»(()  * 

+  A„(«())  +  A'(P(I))^  an(l)  +  -^Xn(t)dt 

(f) 

VL0 

vc0 

+‘^=r-^£^o(^(0)ao(0^^Ao  (0  +  V  (^(0)an  (0  (0 

*  _ _ _  f  i  ^ 

+_l^^e£o(/*(0)ao(0[l  _  ^(0]  (0  +  0  “Tp^ 

Vco  W*-o  y 


(A.3) 
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dv(t)  +  dXc(t)  — 

um)+W))^m  a0(t)+-jLx0(t)  i-V(t)--jLxc(t)dt 
+  «n(W ))+&(p(t))-jLY(t)  an(t)  +  -j=Xn(t)  i-v(t)--jLxc(t)  it 
+^r^e£o{P{t))ao(t)[i-  ~  *K0]  ^  (0 


-  »K0]  dW^  (t)  +  of  4 


(A.4) 


^  =%«(0/«o(0,^(0^(0) 


+5^  (a„  (f  ),a0  (0/  W  T(0) 

+8'an(an(t),a0(t),mMt ))^S 
+$atf(a»(0/M0'/*(0'T(0)^^ 

+-7^=^r(W«n(0/«o(4‘r(0)dw< 7T(0 

VLo 

1  fll 

+-7—- aTrfWT(f)+o  -t— - 
VLo  WLo ; 

The  terms  of  order  1  in  (A.l)  -  (A.5)  yield  the  deterministic  equations  (5.5)  - 


(5.9). 
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The  terms  of  order  l/\/Co  in  (A.l)-  (A.5)  give  the  following  stochastic 
differential  equations  for  the  noise  terms. 

-en(p(t))xg(t)^  -  no 

-S.W0M0 

-^nowo-woHwno 

+M(/5(f)|xc(()-X„(()-X0(()-Xs(l)] 

+/nno)Mo  - «» w  - «» w  -  noino  <a.6) 

-Jrt„(P(l))an(t)  dW^(t)-  ^e0(y3(())«„(f)  <J  W0(f ) 

+^£2ju(/3(t))(T}(t)  -  an  (t)  -  a0(t)  -  y(t))  dWM(t ) 
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= -umKM-KmMW) 

at 

-0.(«OMO 

HW))xn(i)+*\mMWf) 

-je^rnKfi)  *wk  w + -Mmkw 
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d-^p. = -enm)xg(t)^  - 

-mm  x‘« 

-mm 

-tnWQKlW) 

Hn  wtw  -  nm„(t) + m(‘)Km  -  iM)n<) 

-UmKWM + &(«'))(!-  iKO)Xo(‘) 

+«(/3(f)K(,)(i-n(f))ni) 

+A0(/3(l))X0(i)  +  A;(p(f))a0(l)nO 

(A.8) 

+An(/?(0)xM(f)+A;O(f)K(0n0 

-j£0„(P(t))r(t)a„(t)/rj(t)  d\Vdn  {t)  -  Jee0{P{t))Y{t)a0{t)/Ti{t)  dW6o  {t ) 

+^0(P{t))a0{t)  dVJXo  (0  +  V^nW0K(0  dwK  W 
+J^n(/3(0)a„(f)[i-^(0]  <^n  (0 + ^0{p(t))ao(tp-n(t)]  dw^  (t) 
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(A.9) 


= -&(«())«»  (<)xc(<) + £(/3(<)H(<)[i  -  nMM') 
+&(«0)[i-n(0]x„(<) 

-umKWcit) + awoW'X1  -  imi*m 

+£»(/5(0)M(f)]x„(f) 

+1/^„(/!(i))a0(i)[l-T!(()]‘**V4  (t) 
+^„(/J(0)a„(()[l-n(()]‘iW4(t) 

^  =  5^(^((),a„(<),a„((),T(<))y(f) 

+5«o  (^(f),a„(0,«o(0^(0)xo(0 

+CTT(/3(f),an(0,ao(0^(0)^war(0 
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